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ABSTRACT 

The problem of separated flow about a circular cylinder is studied 
analytically through the use of two potential flow models. Following 
a detailed review of the present state of knowledge, mathematical 
formulation of the models and appropriate computer programs are 
presented. These models are shown to be capable of simulating the 
asymmetric vortex shedding. 
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TABLES OF SYMBOLS AND ABBREVIATIONS 



ENGLISH 



C D 


Drag coefficient; defined by equation (1) 


C P 


Pressure coefficient; defined by equation (27) 


c 


Radius of circular cylinder or sphere 


D 


Differential operator 


d t 

f d 


Cartesian tensor notation for Eulerian total derivative 
Drag force on circular cylinder 


f l 


Lift force on circular cylinder 


f . 
1 


Cartesian tensor notation for body forces 


g 


Acceleration of gravity 


i 




K 


A free parameter in Free Streamline Theory 


M 


Mach number or number of free vortices 


N 


Total number of free and attached vortices in a 
flow field 


P 


Pressure 


n . 

i 


Cartesian tensor notation for the three axial velocities; 
q i = u; q 2 = v ’ and ^3 = w 


Re 


2 . c 

Reynolds number; for circular cylinder, Re = ^ 


S 


2 C 

Strouhal number; for circular cylinder, S = 

U T 

OO 


T 


Strouhal period; T = length of wake generated in 

one period 


t 


Time 


U 


Velocity of the ambient flow; could be time dependent 


u 


Instantaneous velocity in x-direction of a two or three 
dimensional cartesian system 
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V 



Instantaneous total velocity vector 



v 



W(z) 



w 



Z 



GREEK 

r 

€. . , 
JL+ A 
«/ 

V 

6 

H' 

Y 

f 

4 

i 1 

SUBSCRIPTS 

i 

1 

k 

m 

n 



Instantaneous velocity in y-direction of a two or 
three dimensional cartesian system 

A function of complex variables; defined by ; 

W(z)= 'l' (x , y) + 4 > (x,y) 

Instantaneous velocity in z-direction of a three 
dimensional cartesian system 

Complex variable; z = x + iy, an undefined location; for 

example, z = x + iy could be the location of the nth 
n n J n 



Strength or circulation of a vortex; defined as the line 
integral of the tangential velocity component around a 
closed curve, or defined as the integral of vorticity 
over a surface bounded by a closed curve 

Specific weight; defined by, $ ~ p $ 

Cartesian tensor permutation symbol 

Cartesian tensor vorticity vector; defined by equation 
(20) in Section II 



An angle measured from the positive x-axis 
Absolute viscosity 



Kinematic viscosity; defined by, r 
Density 



= M- 

t 



Velocity potential function 
Stream function 



Cartes 
quant i 

Cartes 

Cartes 

Vor tex 

Vortex 



ian tensor index; 
ty pertains to an 

ian tensor index; 

ian tensor index; 

counter index; m 

counter index; n 



i = 1,2,3; or indicates 
image vortex 

J = 1,2,3 

k = 1,2,3 

= 1,2, ,M 

= 1,2, N 
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I. INTRODUCTION 



The postulation of an analytic solution for a viscous fluid flow 
about a circular cylinder has perplexed hydrodynamic investigators for 
several centuries. Although notes left by Leonardo da Vinci (1452-1519) 
indicate more than a qualitative interest, it was not until the early 
seventeenth century, during the evolution of Potential Theory, that the 
first formal attempts were made at quantitative analysis. Though there 
was some understanding of viscous effects during this period, they were 
basically qualitative and were not taken into consideration when in 1750 
D* Alembert found that: 



This result, classically known as D'Alembert's Paradox, was of course 
unsatisfactory. Decades have elapsed since investigators found the 
qualitative reasons for the erroneous results given in this and ether 
cases by Potential Theory. 

Philosophically speaking, to the engineer the utlimate goal of any 
analytic procedure is to provide a power of prediction for design ap- 
plications. As yet this capability is generally only realizable from 
experimental investigations and then only for certain flow characteris- 
tics. For a circular cylinder immersed in a steady cross-flow as shown 
in Figure 1, the drag force can be calculated from: 



The force on a body in steady motion relative to an other- 
wise undisturbed ideal fluid is zero. [1]* 




( 1 ) 



1 



The bracketed numbers refer to the item number in the List of 



References . 



9 



IV 




Figure 1 . Circular Cylinder lone reed in a Steady Fluid Cross - 
Flow and Indicating Relationship to a Cony lex 
Coordinate tyeten. 
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Cp is determined experimentally as a function of the Reynolds number as 
shown in Figure 2. Except for Re < 1, there is presently no general 
analytical method available to generate these results. This exception 
will be discussed in more detail in Section II. The primary point in- 
volved here is that for each geometric shape, or group of geometric 
shapes immersed in a steady flow field, there would have to be an ex- 
perimentally found curve for each case. In fact, for an incompres- 
sible fluid, the form of equation (1) would have to be more general: 

F d = f (Re) , (2) 

and for a compressible fluid: 

F d = f (Re , M). (3) 

For all but the simplest geometries, the definition of Re and M are in 
themselves difficult problems; not as to their general meanings, but as 
to specific combination of the physical parameters for a given applica- 
tion. 

Though used only as an illustrative example above, the analytical 
determination of has more subtle implications than those previously 
indicated. It requires a complete understanding of separated flow, 
particularly from the aspect of vorticity generation in the boundary 
layers before separation. It is this point that presents the greatest 
problem in the postulation of a mathematical model; for it is assumed 
that if a general analytical solution can be found for in time- 
dependent flows, then analytic determination of all or most of the other 
flow characteristics will follow. Also, since a circular cylinder and 
its surrounding flow field can be easily represented in two dimensions 
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leurp 2. Lop -lop Flot of Drsp Coefficient Versus Reynolds Number tor 
a Circular Cylinder Irwnersed in s Stssdy Cross-Flow. |2| 



by a function of complex variables, it is possible to easily transform 
the cylinder into other two-dimensional shapes using the technique of 
conformal mapping. For inviscid, incompressible fluids, this technique 
has been used to great advantage particularly in aerodynamics, however 
external knowledge of boundary layer effects was needed. 

So far, the discussion has centered around total drag. For future 
analysis it is necessary to break this into its two classical components, 
skin friction and form drag. The skin friction is the integral of the 
shearing stresses over the surface of the body and is dependent on the 
action of the fluid in the boundary layer. The form drag is equal to 
the integral of the forces normal to the surface of the body and is 
assumed to be dependent only on the inviscid fluid action outside the 
boundary layer. This assumption is made since an order of magnitude 
analysis of the Navier-S tokes equations indicates that the pressure 
gradient through the boundary layer is negligible as compared to other 
effects. The practice of considering the fluid outside the boundary 
layer as inviscid is due again to an order of magnitude analysis of the 
Navier-Stokes equations. This time, when considering only the region 
outside the boundary layer, the equations reduce to the inviscid equa- 
tions of motion and permit the application of Potential Theory to this 
region. Note that large Reynolds numbers have been assumed in both 
instances . 

Following the above reasoning further, leads to the conclusion that 
if no separation of the boundary layer occurs, the form drag is essen- 
tially zero since there is no disturbance to the external inviscid flow. 
At very low Reynolds numbers, this conclusion is reasonable and can be 
considered a basis for the famous "creeping motion" solution of G. G. 
Stokes in 1851. In this case skin friction predominates. At larger 
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Reynolds numbers, however, there is boundary layer separation and form 
drag becomes important. The external flow field is disturbed by regions 
of intense vorticity which causes a decrease in pressure behind the body 
as compared to the undisturbed inviscid solution. It is evident then, 
that without considering the direct shear force effects of viscosity as 
represented by skin friction, and the indirect effects represented by 
vorticity regions, the qualitative concepts reduce such that D'Alembert's 
quantitative solution results. 

Earlier in this discussion, mention was made concerning the practice 
of considering the undisturbed fluid outside the boundary layer as in- 
viscid. This requires amplification and a word of caution. In 1908, 

H. Blasius presented a solution to the boundary layer on a circular 
cylinder. As far as his paper was concerned the method of solution, 
rather than the result, was the most important contribution. Briefly, 

it involves the use of a series expansion for the free stream velocity 

2 

just external to the boundary layer. If this velocity is taken as that 

found for incompressible, inviscid flows derived from Potential Theory 

on the surface of the cylinder, the angular distance between the upstream 

stagnation point, and the boundary layer separation point is found to be 

3 

108.8 degrees; the center of the cylinder is taken as the origin. If 
experimental velocity determinations are used, the solution yields a 
separation point at about 82 degrees. The observed separation point is 



For a complete presentation of this method, see [2], pages 154-162. 

The boundary layer separation point is defined as: ^u-c 

where u is the fluid velocity parallel to the body surface, 
and y is directed normal to the surface, with its origin at the surface. 
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at about 81 degrees. Herein then lies the problem with the above prac- 
tice: even at high Reynolds numbers there is sufficient interaction 

between the boundary layer and the external "inviscid" flow to cause 
considerable error in the results of essentially exact methods of solu- 
tion. As will be shown in Section III, the separation point is critical 
in any method that attempts to mathematically generate a reasonably cor- 
rect drag coefficient curve. 

In recent years, extensive experimental work with the flow past 

bluff bodies has provided researchers with a clearer understanding of 

the mechanism of boundary layer separation. As a result, several models 

have been devised to portray the overall fluid action. However, to date 

4 

none has provided the sought-after general solution. A survey of the 
literature indicates that there are two basically different analytical 
approaches being made. The first, and probably the most obvious, is the 
solution of the nonlinear Navier-Stokes equations using finite differ- 
ence techniques. The second is a potential flow model in which the shear 
layers are represented as point concentrated vorticies. Both these 
models, particularly the latter, will be described in greater detail in 
succeeding sections. 

Since this thesis does not propose a general solution, it is hoped 
that the following sections will help to clarify the problem and to 
provide a basis from which future work can be done. Section II reviews 
the literature on separated flows; concentrating on those aspects that 
relate to incompressible flow about circular cylinders. Section III 
presents the two potential flow models that the author used in an at- 
tempt to generate the aforementioned drag curves. 



The development of the high speed, large storage capacity digital 
computer has also been extremely important to the study. 
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II. REVIEW OF THE PRESENT STATE OF KNOWLEDGE 



Because of the large number of references available on separated 
flows, the author has elected to limit this review through the use of 
the following guidelines. First, since there still exists a lack of 
basic understanding of the mechanism of separation only references to 
studies involving essentially incompressible fluid will be used. This 
enables the added complications of compressible phenomenon, both qualita- 
tive and mathematical, to be ignored. Secondly, a study of the literature 
reveals that several different geometries have been used both in experi- 
mental and analytical work; the results not always consistent nor related 
to each other. Therefore, the circular cylinder is taken as the basic 
geometry for this review. The use of a circular cylinder implies, to a 
large extent, two-dimensionality; which is another simplification. 

Though particularly helpful from a qualitative standpoint, this simplifi- 
cation is also useful mathematically. Thirdly, the information will not 
necessarily be presented in chronological order as is so often done in 
reviews of this nature. It is felt that a loose division by Reynolds 
number is more pertinent since the analytical solutions thus far derived 
are only applicable to given ranges of this M magic number," as von 
Karman calls it. Lastly, the discussion will in general be limited to 
the laminar flow regime. As the Reynolds number increases above the 
critical value, increasing consideration must be given to compressibility 
and eventually to cavitation; two complicating subjects outside the scope 
of this thesis. 

A . REYNOLDS NUMBER LESS THAN ONE 

In this Reynolds number range, the viscous forces acting on the 
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cylinder are considerably larger than the inertial forces. This leads 
to the conclusion that if the inertia terms are dropped from the Navier- 



Stokes equations, an approximate solution may result that would compare 
favorably with experimental data. Stokes, as previously indicated, did 
this in 1851 when he presented the first solution for the drag on any 
body. Though he used a sphere in his solution, it is still worthwhile 
presenting here. The mathematical details are not the same, but the 
basic steps are. 

The complete Navier-Stokes equations for incompressible, constant 
temperature conditions are (expressed in three-dimensional Cartesian 
tensor notation [4]): 



By dropping the inertia and body force terms, then cancelling density 
from both sides of the rearranged equation, one has 



The four equations represented by (5) and (6) have the same usual 
boundary conditions used for any boundary layer problem of this nature 
when the body is immersed in an infinite fluid: the fluid velocities 

normal to and tangent to the body surface are equal to zero. 

After taking the divergence of both sides of equation (6), it 
becomes : 




The incompressible continuity equation is: 




(5) 




( 6 ) 




( 7 ) 
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But the right hand side of equation (7) can be written as 



A (ft, L ) 



))) 



which is the second derivative of the continuity equation and there- 
fore equal to zero. This means that in creeping motion, the pressure 
field satisfies Laplace's equation and hence is a potential or harmonic 
function. Now, by introducing the two-dimensional stream function. 








1 



into the right hand side of equation (7), dropping 
two-dimensional biharmonic equation results. 




( 8 ) 

and JX the 



D u ^ = D 2 (J > 2 f) 



Without going into the mathematical details, it is sufficient to 
say that the pressure and velocity components reduce to: 



u = % 




V = 






vV = 




C *2 

r 3 





2 ? * 2 

r - x +- f z 



( 10 ) 



The geometric configuration is such that the center of the sphere is at 
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the origin of the coordinate system and is placed in a positive velocity 

field, u , which is parallel to the x-axis. 

oo 

Next, by using the standard definitions of shear stress along with 
equation (10), the viscous and pressure forces can be found by integra- 
tion over the sphere 1 s surface. The total drag for the sphere then 
becomes , 



Fj) = £ ^ M c &oo ' 



(id 



If a drag coefficient is formed in the normal manner (i.e., referring 
to the dynamic head and frontal area), then. 




( 12 ) 



Equating (11) and (12), then, 

Z± 

Re 




(13) 



Comparison of this result with experimental data indicates that 
these equations are valid only for cases where Re K 1. A classic case 
of application is that of small, light particles, assumed spherical in 
shape, settling in a fluid. Assuming that the particle is at its ter- 
minal velocity, then summing the forces yields. 






(14) 



Letting A r- tijO sph 



ere 



which is the famous Stokes 



fluid^ 
1 number. 



then equation (14) reduces to 




(15) 
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Kieure i. Flow Line Formation Around a Circular Cylinder for 
3 < Re < 40. Bottom is mirror image of this. 




C. W. Oseen, in a paper published in 1910, improved on the above by 
taking the inertia terms of the Navier-Stokes equations partly into ac- 
count [2]. By doing so, the resultant drag coefficient became: 



Comparison of (16) with experimental data shows favorable correlation 
up to a Reynolds number of about 5. 

In recent years, other investigators have formulated other solutions, 
however, their main contributions have not necessarily been to engine- 
ering but more to mathematical prowess. Kaplun [5] based his on an 
asymptotic solution of a series expansion of the Navier-Stokes equations. 
The results are generally restricted to two-dimensional bodies with 
equivalent Reynolds numbers less than one. Sir Horace Lamb [6] also 
presented a solution applicable to this range. He followed the basic 
ideas of Oseen, however, where Oseen used a sphere, Lamb used a circular 
cylinder. 

B. REYNOLDS NUMBER BETWEEN THREE AND FORTY 

In this Reynolds number region, the boundary layer is still laminar, 
however, separation occurs causing the formation of two vortex cores in the 
wake of the cylinder as indicated in Figure 3. This region is sometimes 
referred to as the "Twin-Vortex Stage. 1 ' A great deal of experimental evi- 
dence indicates that the cores are generally stable in this region pro- 
vided that no external excitation occurs. It is in this region also 
that the greatest successes have occurred with the "exact" solution of 
the Navier-Stokes equations. Remembering that, above a Reynolds number of 
about one, the inertia forces are too large to neglect, the phenomenon 




# 



(16) 



21 



This is most 



requires the solution of the full Navier-Stokes equations, 
readily done through the use of either the forward or backward finite 
difference techniques. Until recently, without the advent of the digital 
computer, these methods have been difficult to employ. Thom [7] in 1933 
used a relaxation method of hand calculation to produce the streamline 
pattern shown in Figure 4 at a Reynolds number of 20. 

It is true that the relaxation method can be programmed for a digi- 
tal computer and made to converge to any realistic degree of accuracy, 
however, it is also true that computation time increases rapidly for a 
slight increase in desired accuracy. Thus this method is not considered 
the most desirable. The Gauss-Siedel iteration method is somewhat 
different in procedure, but can require the same amount of computer 
time for solution. It should be noted that both of these methods can 
be modified if done by hand since one can make shortcuts that are diffi- 
cult to program for the computer. Even so, the digital computer reduces 
the time required from " about one year and a half with twenty work- 

ing hours every week, with a considerable amount of labour and endurance 

"[8] to a total computation time of about three to five hours. 

Even this could probably be reduced by streamlining the iterative proce- 
dures. Robertson [1, p. 366] summarizes these methods quite well. 

Though these analytical methods produce interesting flow patterns, 
their primary reason for use is to give needed information about sepa- 
rated flow. Unfortunately, they do not appear to explain the mechanism 
of separation or more importantly the feedback aspects of the wake on 



Body forces, f . , are usually neglected in these solutions since 
their presence adds nothing to the understanding of the problem and in 
general add unwanted complications to the mathematics. 
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at He * 20, by Thom (1933) using Relaxation Method. 
The Undisturbed Flow Streamlines are Represented 




the vortex generation in the boundary layer. The experimental and ana- 
lytical work does yield approximate Reynolds numbers for the beginnings 
of separation and the formation of the Karman vortex street; about 5 and 
40 respectively. From an analytical standpoint, the mesh size of the 
finite difference grid and the method of defining the conditions along 
the outside grid points^ can affect the resultant separation and wake. 
Experimentally, just the placement of the confining walls of the ap- 
paratus used for testing can increase the Reynolds number, at which the 
wake becomes unstable, by more than 300 percent [9]. It should be noted 
here that although the wake is generally considered to be stable for 
steady flows in the Reynolds number range under discussion, it can be 
made unstable both theoretically and experimentally under artificial 
excitation. The ease with which this can be done, of course, increases 
as Re approaches the range of instability, making exact determination 
of the Reynolds number at which the vortex shedding begins very diffi- 
cult. This same idea holds true for the determination of the separation 
Reynolds number. In both cases there appears to be some experimental 
evidence to indicate that the Reynolds number is not the exact similarity 
parameter that should be used in comparing results; the reasons for this 
are not entirely known. It is the belief of the author, however, that 
these disparities in results are probably due to variance in technique 
rather than to any fault of the parameter. It is recognized that in 
actuality there may be a three-dimensional effect that has not been taken 
into consideration theoretically or experimentally. Herein may be much of 
the reason for the differences in experimental results. 



Requiring that the velocity on the surface of the cylinder be 
zero provides the boundary condition for the grid points on the 
cylinder surface. 
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Figure 5. Beginning* of Vortex Street, Re ^ 40. 




Figure 6. Kama* Vertex Street Indicating Cyclic 
•Hedding Vortex Cores, 80 < Re < 300. 
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C. REYNOLDS NUMBER BETWEEN FORTY AND THREE HUNDRED 

The primary phenomenon represented in this region is the Karman 
vortex street; from the condition of an almost stable wake with slow 
viscous decay to the ultimate condition of cyclic detachment of vortex 
cores. Figures 5 and 6 indicate the flow patterns at these two "ex- 
tremes." This region is also characterized by a smooth decrease in drag 
coefficient and an almost linear increase of Strouhal number, S. 

Though study in this region has not produced any startling insight 
into the mechanism of separation, it has provided much information about 
vortex streets. Basically, it has been found that given two parallel 
vortex sheets separated by some initial distance, a solid body need not 
be present for a vortex street to develop [10]. The creation of these 
sheets is of course dependent on the viscous boundary layer effects 
caused by the presence of the solid body. 

Due to the oscillatory nature of the vortex core separation into 
the wake, there is a cyclic shift in the separation points of the 
boundary layers on either side of the cylinder. This action is evidence 
of the strong "feedback" effect that the wake conditions have on the 
boundary layer in this Reynolds number region. This condition exists 
during the so-called steady state period after the various parameters 
such as drag coefficient, Strouhal number, etc., have settled to fairly 
constant values. In impulsively started flow situations, the flow pat- 
terns above and below the cylinder are basically mirror images of each 
other. The mechanism that causes eventual asymmetry is not clearly 
understood. In recent years, with the help of the digital computer, 
this asymmetry has been generated in finite difference solutions by the 
artificial insemination of small amounts of extra vorticity into one of 
the boundary layer separation regions, once the twin core pattern was 
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fairlv well developed. Using this method, Fromm and Harlow [11] were 
able to duplicate the flow patterns around rectangular cylinders at 
these intermediate Reynolds numbers. At the lower numbers the system 
"absorbed" the extra vorticity with only minor disturbances to the wake, 
while at higher Reynolds numbers, the vortex cores separated cyclicly to 
form a vortex street. 

Though the details of solution of the Navier-Stokes equations varies 
with the investigator, the same general procedures are used by all. The 
following is an attempt to clarify the salient points involved, but 
first, the reason for using these equations should be reiterated. 

Through this Reynolds number region, there appears to be a balance of 
sorts between the viscous and inertia forces, or effects, which implies 
that neither can be neglected or ignored. In the potential model that 
will be discussed later, the viscous forces are neglected from a drag 
standpoint. In the Stokes solution for Re ^ 1, the inertia forces were 
neglected with no apparent loss of practical accuracy. This region then 
appears to be nature 1 s transition region from one extreme to the next, 
and therefore the complete governing equations must be used. 

The method, in general, is as follows; dropping the body force 
term from equation (4), and expanding to the two-dimensional case, 
there results: 




(17) 



and , 




(18) 
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Eliminating pressure by cross -dif ferentiation and subtraction, one has: 
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The vorticity vector is defined as 






which means that since only two dimensions are involved, vorticity 
exists only parallel to the Z-axis. Therefore, letting 



( 20 ) 



y_ y _ _ yu 

^ ^3 ~ 



(21) 



and substituting back into equation (19), one has 



, Y dA 4. Y 4 2^ — + V ~ - / ’ [ - ^ ~ 4 ~ ^ ^ 1 

w + + u *x v a* ~ ' + w-l 









or finally. 



a % 4 . ^ av S 
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c) 
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- T Lw\/ 



( 22 ) 



u, v, and P must be assumed continuous functions of x, y, and t. 
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The system of equations needed becomes complete by introducing 
the stream function relationship. 



^'3 



(8) Repeat 



into the continuity equation, = For t * ie i rrotat i° na l case, 

this would result with the Laplacian of the stream function equal to 
zero. However, there is rotation present and this introduction must be 
handled differently. By adding ^ to both sides of the continuity 
equation, one obtains 



+ iy _ Zu _ y , 

c)x o)X c)X 3 

The rearranging of the left side produces; 



(23) 



[u-hvl + i- [v -ul = * < 24 > 

c)X L L ° 



Now substituting for u and v the appropriate expanded relation from 



equation (8), i.e.. 




and v = 






, and then performing 



the partial differentiation as indicated in front of each bracket, the 
equation becomes: 



_ tl- = % 

dx z Jxdy. dj / 2 



This of course reduces to a form of Poisson's equation, 

XV XV _ _r 

Jx 2 if ~ * ' 



( 25 ) 



Equation (22), when written in terms of ^ and ^ reduces to 



which is somewhat easier to program since the necessity of finding the 
velocity gradients has been eliminated. If desired, these could be 
found from the other forms of the equations given above. As before, 
the assumption concerning continuous functions noted in footnote 7 still 
applies . 

The specific finite difference formulas used to convert equations 
(25) and (26) into a numerically solvable system depends basically on 
the largest allowable truncation error and on the mesh size of the grid 
points, Crandall [12] discusses quite thoroughly all the aspects in- 
volved for both these considerations. It is sufficient to say that 
greater accuracy requires less truncation error and a finer mesh size. 
Generally, since finer mesh sizes require larger computer storage cap- 
abilities, they are used only in the areas of the flow field where large 
gradients are expected. 

There are in general two basic methods of solving this system of 
equations. The one most widely used, until recently, has been of an 
iterative nature such as the relaxation or Gauss-Siedel methods; or 
modifications thereto. The primary concern here is that of convergence 
and stability of solution. The Gauss-Siedel method is adversely sus- 
ceptible to both of these while the relaxation method generally does 
not have a stability problem, but can require an abnormally large number 
of iterations before some degree of convergence is attained. These 
methods are however, fairly simple to program. With the advent of the 
digital computer, more sophisticated numerical methods can be used to 
solve the system of equations explicitly using matrices. Methods such 
as Gaussian elimination with back substitution can produce a very ac- 
curate inverse for a matrix that is not singular, and a properly 
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constituted matrix representing a physical system should not be singular. 
For large systems of equations, the coefficient matrix can be reduced to a 
triangular form through the use of a sequence of linear operations. The 
inverse of this reduced matrix is then used rather than the inverse of 
the complete system. 

Besides Fromm and Harlow [11] several other investigators have 
studied this Reynolds number region using modifications of the above 
procedures. Most numerical solutions have involved rectangular cylinders 
due to the ease with which Cartesian coordinate finite difference equa- 
tions can be applied. Thompson [13] used a Gauss-Seidel method but 
appears to have assured convergence through the use of an "acce leration 
parameter" which is determined from a consideration of the eigenvalues 
of the matrices involved. His results compare favorably with experi- 
mental data, however, some of his computer runs required as much as nine 
hours on an IBM 360/40. This could have serious implications as far as 
accuracy is concerned, since each time step must adequately converge be- 
fore the conditions for the next time step can be calculated. If con- 
vergence presents no problem, then there is another important result of 
his work that should be mentioned: no artificial introduction of vorti- 

cltv is reauired to cause the eventual formation of the vortex street. 

Wang [14] used cylindrical polar coordinates in a matched asympto- 
tic expansion solution of the Navier-Stokes equations. Though the 
results offer no new insight, the method is indicative of the mathe- 
matical attempts being made. Ingham's paper [15] offers constructive 
criticism of Payne’s method of finite difference solution [16], and is 
indicative of the pitfalls waiting for the unwary. Ingham's examples 
of variation of calculated for different mesh sizes at the same 
Reynolds numer is quite dramatic. 
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Figure ft. Experimental and Caaiputed Pressure 
Coefficient and Computed Surface 
Vorticity (5 ) at Lorn Reynolds Numbers. 

A Computed C p , Apelt, Re - 40 
B Ccaputed Cp, Kavaguti, Re ■ 40 
C Co m puted 5 , Apelt, Re - 40 
D Experimental, Linke, Re * 2800 
E Computed Cp, Potential Flow Theory 
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The pertinent experimental and analytical data derived by investi- 
gators in this Reynolds number range are shown in Figures 7 and 8. 

Important contributors of this information include Taneda [17], Tritton 
[18], Roshko [19], Rosenhead [20], Wille [21], and Apelt [22]; besides 
those previously mentioned. 

D. REYNOLDS NUMBER BETWEEN THREE HUNDRED AND CRITICAL 

This subcritical region is in general characterized by a relatively 
constant mean drag coefficient and a progressive loss in the distinct 
visualization of the Karman street. As the Reynolds number increases 
there appears to be less "feedback" effect of the wake on the separating 
boundary layers which are still laminar. The Strouhal number has a 
tendency to decrease slightly with an increase of Re. It is in this 
region that the potential flow point vortex model of the feeding sheet 
appears to have the best possibility of providing information since the 
inertia forces are so much larger than the viscous forces. The use of 
the finite difference Navie r-Stokes equations to provide a solution 
seems to break down due to the increased turbulence in the wake. Al- 
though a much finer mesh size might alleviate this problem, the fact 
that the equations now being used do not provide for velocity perturba- 
tion is probably the more basic reason. There appears to be no argu- 
ment that the equations are applicable for use in describing the boundary 
layer on the front part of the cylinder almost up to the critical Re of 
about 5 x 10^. As the Reynolds number increases from the low end of 
this region, the boundary layer effects appear to feed vorticity into 
the wake at an ever increasing rate causing the wake to become extremely 
turbulent. As this occurs, the well ordered vortex core regions of the 
Karman street essentially disappear. Although the Strouhal number can 
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still be measured, the dominant periodicity seen in the lower Reynolds 
number wake flows all but disappears. Experimental studies using turbu- 
lent flow measuring techniques have confirmed the increasing turbulence 
with increasing Reynolds number. Ironically, the feedback problem di- 
minishes, but the problem of transition from laminar to turbulent flow 
conditions becomes more important. The understanding of the latter 
problem is even less than that of the former. The end result is that 
the applicability of the laminar Navier-Stokes equations decreases with 
increasing Re. 

In the "critical 11 Reynolds number region, loosely defined as 2 x 10"* 
< Re < 10 6 , there is a sharp decrease in the mean down to about 

g 

0.3. According to Schlichting [2, p. 39], this is due to the following 
reason: when the total drag on any body is primarily form drag, the 

transition from a laminar to a turbulent boundary layer causes a decrease 
in the size of the wake and results in a decrease in drag. This can be 
interpreted as saying that the drop in is a boundary layer phenomenon, 
though it does not quite explain the entire mechanism. After this transi- 
tion occurs, increases to a value slightly lower than that before 
transition. There is also a return of a distinct vortex street. This 
time however, the "vortex cores" are highly turbulent internally. As 
the Reynolds number increases beyond this region, compressibility ef- 
fects become greater, further complicating the problem. Depending on 
the fluid, very high Reynolds numbers can bring on the added, sometimes 
predominant, problem of cavitation. 



This is true for a hydraulically smooth cylinder. The roughness 
of the cylinder surface can have a major effect on this minimum mean 
drag and the Reynolds number at which transition occurs. See [2, 

p. 621]. 
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Q % De/rees measured from forward stagnation point* 



Figure Q . Experimental Pressure Coefficients Found 
at Medium and High Reynolds Numbers 
Indicating Laminar (S^) and Turbulent 
(Sj) Separation Points. 

A Linke, Re » 2.8x10-* D Flachabart, Re * 6.7x10^ 



BjVage and? Re * 1.1x10^ 
CfFalknerJ Re - 2.1xl0 5 



X Roehko, Re • 8.4x10^ 

F Potential Flow Theory 
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Experimental studies of the pressure distribution on the cylinder's 
surface have yielded some valuable information necessary for future 
analytical postulation. The pressure coefficient, defined as 



r ~ % AP 

S 3 yO u z 

/ oo 



(27) 



is an excellent parameter for theoretical and experimental comparison 
particularly for creeping motions (Re < 1), and at Reynolds numbers where 
the boundary layer separation points maintain fairly constant positions. 

At very low Reynolds numbers, the comparison is quite good especially on 
forward part of the cylinder, between points A and B in Figure 1. Be- 
tween points B and C, the actual pressure decreases somewhat since energy 
is expended in the boundary layer; i.e., not all the kinetic energy 
acquired from A to B is recovered as potential energy as the fluid de- 
celerates toward C. In Potential Flow Theory, points A and C are of 
course stagnation points, and whereas the pressure at A can be well ap- 
proximated by this theory over the entire Reynolds number range, the 
pressure at point C cannot as indicated in Figure 9. Note how curve D, 
found for a Reynolds number above the critical, more nearly approxi- 
mates curve F (potential flow). Photographs of the wake pattern pro- 
duced for this post-critical Reynolds number shows a narrow wake with 
the rest of the flow field similar in appearance to that of an ideal 
fluid. Note also that after the boundary layer has separated, as- 
sumes a nearly constant value, not the same for all Reynolds numbers 
however. 

Though information about the Strouhal number is still incomplete 
above Re 10^ , there is an interesting similarity between the inverse 

of this parameter and C^. The curves for the two values, plotted over 
the whole range of Re from the beginning of core shedding through the 
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critical region, are very similar in shape. In fact, the similarity is 
so good that the drag coefficient could be written as 




with the resulting value having an accuracy within about 15 percent. 

The explanation for this result is not complete, however, the point here 
is that there is apparently a strong conceptual connection between core 
shedding frequency and drag. It is interesting to also note that through 
the region of nearly constant drag, not only is 1/S fairly constant, but 
that the oscillations of the boundary layer separation points are gen- 
erally quite small in amplitude. 

Investigators have attempted several approaches in the search for 
an analytic solution in this region and there have been some notable 
successes. Each approach, however, involves the use of at least one 
experimentally determined or arbitrarily selected parameter. 

The free-streamline theory developed by Kirchoff, idealizes the 
free shear layers that separate from the cylinder by considering them to 
be surfaces of velocity discontinuity; hence a free streamline. This 
supposition works quite well for liquid jets discharging into a gas, 
however, for cylinder wakes, because of the assumption that the pressure 
in the wake is the same as the free stream pressure, the resultant cal- 
culated drag coefficients are much lower than those experimentally observed. 
The agreement can be greatly enhanced by the introduction of an arbitrar- 
ily selected parameter usually symbolized by K. This parameter is con- 
cerned with the depth of discontinuity (notch) in the logarithmic 
9 

hodograph plane, and is therefore concerned with the ratio between the 



In this plane, the two-dimensional physical flow is represented by 
straight lines parallel to the two coordinate axes, and is the result of 
conformal mapping. See [1, p. 351]. 
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free stream-line velocity and the free stream velocity. For very high 
Reynolds number flows, where cavitation exists, K is fixed by the cavi- 
tation index, CT . Fairly good agreement between experiment and theory 
exists in this case, particularly with respect to cavity flows. 

Dennis and Shimshoni [26] have probably been two of the more suc- 
cessful investigators in analytically deriving the mean drag coefficient 
over the entire laminar boundary layer Reynolds number range. Their solu- 
tion involves series expansions^ of the Navier-Stokes equations, set up 
for finite difference evaluations. At slow, creeping motions, their 
solution converges to that of Stokes and Lamb, but at Reynolds numbers 
above about 1.5, the values of are about 5 to 20 percent higher. 

They allow no asymmetry of vortex shedding. The mesh size they use ap- 
pears to be much smaller than that of other attempts. 

As noted earlier in the discussion, this Reynolds number region is 
characterized by the predominance of inertia forces over viscous forces. 
Also noted was that this leads to the possibility of using point con- 
centrated, potential flow vortices to represent the shear layers or 
feeding sheets. Among those who have worked with this approach, 

Sarpkaya [35], Bryson [29], and Gerrard [25] have made some notable 
strides. The description of the work of the first two will be left to 
Section III, since it is basically with their ideas that the author 
worked. Gerrard* s work, though similiar in mathematical approach to 
the others, appears to indicate a greater concern for wake characteris- 
tics than for finding solutions to the basic separation problem. 
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Similar to Oseen and Lamb. 
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Besides the basic points covered here, there is a great wealth 
of experimental information on the more detailed aspects of the sub- 
ject that should be covered by any reader desiring to proceed further 
with the investigation. Toward this end, the author has attempted to 
compile the more notable of these in the List of References. Although 
some references have not been mentioned specifically herein thus far, 
this is not intended to mean they are of lesser value. 
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III. TWO ANALYTIC MODELS 



A. INTRODUCTION 

In the previous two sections, the discussion has centered on those 
few basic known facts and some of the many problems involved in this 
study. This section will present, in detail, two different analytic 
models, based on Potential Flow Theory, which are theoretically ap- 
plicable to the relatively higher Reynolds region discussed in Section 
II (D). 

The first model, called Model A hereafter, is conceptually more en- 
lightening than the second; hence the reason for its presentation first. 
It is basically derived from Bryson's [29] model, but where his model 
was limited to the growth of axisymmetric vortex cores, the equations for 
Model A are recast to allow "forced" asymmetry to be included into the 
solution. To the author's knowledge, this is the first time this has 
been done. This method, like almost all the others, has no provision in- 
cluded in the governing equations for "natural" asymmetry. Thus if asym- 
metry is to be studied, its beginnings must be forced. In brief terms 
then, the basis of this model is that large vortex cores are fed vorti- 
city through an infinitely narrow feeding sheet, or shear layer, con- 
necting the center of the core with assumed stagnation points on the 
cylinder ' s surface . 

In the second model. Model B, the shear layers are represented by 
a string of small strength vortices, called a vortex sheet; see [1, p. 
111]. The equations for this model were also derived and programmed to 
allow forced asymmetry. Whereas Model A represents a continuously fed 
large vortex core, this one creates numerous small vortices which are 
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allowed to roll-up, coalesce and/or cancel at will. There are four birth 
areas in this model, two for the primary vortices, and two for the sec- 
ondary vortices. ^ This arrangement is a modification of Sarpkaya ' s 
model [35] in which there are only two axisymmetric birth areas. Com- 
parison between the two will be left until later in this section. 

B. BASIC POTENTIAL FLOW EQUATIONS 

The basic general equations used in both models will be presented 
in this subsection. The actual application will be left to the dis- 
cussion of that particular model. It is not felt that the history be- 
hind these equations is relevant to this thesis since the interested 
reader can refer to Sarpkaya [36] for this information. 

As mentioned earlier, in Section I, a two-dimensional flow system 
can be represented mathematically by equations written in complex vari- 
ables. Since the mathematical development of this representation is 
quite lengthy, the following is only a brief part of what can be found 
in many reference books and papers on hydrodynamics. To start with, 
it is well to review the basic assumptions. First, the fluid is con- 
sidered inviscid. This implies that there are no viscous forces pre- 
sent and therefore that the models derived should be applicable to the 
higher laminar Reynolds number range. Secondly, since this aforemen- 
tioned Reynolds number range is generally devoid of compressibility ef- 
fects, will be assumed constant. The usual "steady flow n assumption 



Vorticity is shed from the cylinder in the regions of boundary 
Layer separation. These regions are associated with the "primary" 
boundary layers on the cylinder* s upstream side and the "secondary " 
layers due to fluid back-flow on the downstream side. Since vorticity 
is created in these regions, they are called "birth" areas, and the 
small strength vortices created are called "nascent*' vortices. 
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will not be made, because all the flow characteristics appear to be 
dependent on the manner in which the boundary layer generates vorticity 
and on the manner in which it diffuses into the wake. The key problem 
then is the formulation of a method by which vorticity can be generated 
in a supposedly inviscid fluid. This is in direct violation of Kelvin's 
theorem for inviscid fluids which says that vorticity once created, can- 
not be destroyed, and if not already existing, cannot be created. In 
order to overcome this, boundary layer information must be taken into 
consideration. 

From basic considerations of hydrodynamics, an undefined number of 
vortices outside the boundaries of a circular cylinder, suspended in a 
uniform flow field as in Figure 10, can be represented by 




(28) 



There are no symmetry restrictions on this equation. The normally ex- 

12 

eluded terms have been left out. Next it can be shown that for a 
general location Z , on or outside the cylinder, the fluid velocity is 
given by: 



dW(z) 

dz 



-lL-bsiV 




JL 




(29) 



These are the terras that are excluded by virtue of the condition 
that the vortices shed from the cylinder leave a circulation opposite 
to their own on the cylinder. See [38, p. 364]. 
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Fi enre 10. Circular Cylinder in the Pretence 
of Numerous Vortices and a Uniform 
Flow Field. 
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Now that if Z=Z , the velocity goes to infinity. This is to be ex- 
n 

pected since the velocity at the center of an ideal, inviscid vortex 
is equal to infinity. This same equation can, however, be used to 
find the absolute velocity of the vortex with respect to cylinder, here- 
after called the convective vortex velocity, by letting Z®Z and dis- 

n 

regarding the term that goes to infinity. This is due to the fact that 
the motion of a vortex is governed not by its own internal velocity dis- 
tribution, but by the forces induced upon it by other entities in the 
flow field. See [1, p. 124]. 

Next, the drag and lift forces on a body can be calculated from the 
generalized Blasius r therorem. 
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Writing 
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(30) 
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equation (30) can be broken down into its so called "steady state" 
part. 






(32) 



and the time dependent part, 



X z +d.y 3 ~z/>~ jw/z)d 



2 . • 



(33) 



The steady part, using equation (29), can be evaluated [35] such that: 

N 

(34) 
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The evaluation of equation (33) must wait until the discussion of the 
individual models since it involves some assumptions as far as the 
method by which /""* is to vary over the feeding zones. To continue, 
non-dimensionalizing can easily be done by referring all the length, 
velocity and time quantities to and c; or for example, 

Z = cZ' > 



2 1 V — LI ('LL+'tV') ) 

of L 

r = u^r " , 

and , 



Substituting equations (33) into (29); letting U 
and dropping the primes for the sake of simplicity, 



(35) 



= 1.0; c = 1.0, 
this becomes; 
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(36) 



~n 



Up to this point a discussion of image vortices has not been 
necessary, however, certain relations will be needed for future de- 
velopment of equation (33) in succeeding subsections. It can be shown 
bv interpretation of Mi Ine -Thompson 1 s circle theorem [38, p. 157], or 
by using the method of images [38, p. 220], that the second summation 
in equation (29) and consequently in (36), are in reality "image" 
vortices inside the cylinder. The location of the image vortex. 



corresponding to the real vortex at Z is /r7 

n / ~ n 



Its circulation 

is in the opposite sense as compared to that of its corresponding real 
vortex. Positive numerical values for circulation indicate rotation in 
the counterclockwise direction, while negative numerical values indi- 
cate clockwise rotation. In the derivation of the equations, the 



46 



4J <U 



co id 

- >1 

c Id 
O M 

&§■ 

CQ W • 
CO 

<«£ 

O P H 
•H 3 
C 3 01 

0 <u 

01 C M 
H O 

M *H «— I 
id jj id 

a p p 
E »H C 

o o a» 
owe 



cm a 
“Id x 
SOD 
G id P 
O 

CD p/'J-T 

>, Wo o 

M d HH 
00 W X X 



I ID r 

r* 





47 



io asbb-i lo 



symbol has no inherent sign and can therefore take on positive or 



negative values as the situation warrants. 

Since there are conditions well described in [38] that make equa- 
tion (36) not applicable to the internal area of the cylinder, Sarpkaya 
[36] shows that the velocity of the image vortices are strictly a func- 
tion of the position of the image vortex and the corresponding real 
vortex's velocity. By manipulating his expressions, the following equa- 
tion for the image velocity results: 



To reiterate 
when solving 



71. + i V. - 

+ 'L * V 






VL 






(37) 



the term whose denominator equals zero is disregarded 



for - u + iv . 
n n 



C. MODEL A 

As previously mentioned, this model is primarily a modification to 
Bryson's solution [29]. Though his solution appears to give reason- 
able drag coefficients during the beginning time periods of an impul- 
sively started flow condition, once the peak is reached, the drag 
coefficient drops off to fluctuate about zero; Figure 11. This of course 
is not correct as Sarpkaya* s superimposed experimental data [33] indi- 
cates. Though Bryson stops his solution when the vortices begin feeding 

13 

circulation back to the cylinder, this in itself is not enough to 
produce a reasonable drag curve. It was hoped that by studying his 
method, modified to allow asymmetry, more could be learned about the 



The condition occurs when (to be developed later) for a 

given vortex changes sign. If the solution is not stopped at this 
point a negative drag coefficient results as shown in Figure 11. 
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effect of vortex core shedding on the drag coefficient. Before a 
discussion of those results in undertaken, the equations must be 
developed . 

Using Figure 12 as the basis for nomenclature definition, the 
equation development starts with the necessary modifications to equa- 
tion (36). Though these equations can be constructed to allow more 
than two vortices, this will not be done at this time in order to show 
a simplified version of the technique. As will be indicated later, the 
full equations are obtained by the simple addition of a summation term 

to each of the following equations. So, based on two vortices, the 

14 

velocity of vortex 1, if it were free, is given by 



-U +-iv --L -l-h-i-f-S- S. UL 

and likewise for vortex 2; 
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(39) 



-2 Z. 






Next, the feeding aspect of the vortex cores is handled by consider- 
ing the vortex sheet, or feeding layer, as the combination of a single 
concentrated vortex, situated at the center of vorticity, with a con- 
necting sheet to the body of vanishingly small vorticity. This sheet 
can now be considered as a mathematical branch line between Z^ and Z^ 
for vortex 1, and Z^ and Z^ for vortex 2. Since there is a discontinu- 
ity, conceptually, in ^ across these branch lines, there is a velocity 
difference and therefore a pressure difference. This pressure difference 
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A 



free vortex is one that is no longer being fed circulation. 
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Fi tjure 1 2 , Model A, Indicating Capability for 
Asymmetry. 
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Since the dis- 



can be shown (143] and [44]) to be equal to 
tances over which the pressure differences act are ^Z^-Z^Jand^Z^-zJ) 
for vortices 1 and 2 respectively, the forces acting perpendicular to 
these sheets are given by and (Z^-Z^). The 

forces acting directly on the vortices is given, for vortex 1, by 



♦ *,) 



and for vortex 2, by 



. /->/ rfWrf'.Z.) I . i \ 



2 = 2 ^ 

Thus to make the net forces on each sheet and its associated concen- 
trated vortex vanish, ^ the sum of the respective forces must equal 
zero. Therefore, algebraically summing each pressure force with its 
proper vortex force, cancelling out i and , and rearranging the terms, 
the two expressions formed become; for vortex 1: 



Z+/Z-Z , 

' '■ • t> r - d2 _ > 

' 2-Z, 



(40) 



and for vortex 2; 



• <*■> 

These equations indicate that unless is zero, the nth vortex does 

not move with the local fluid velocity, u + iv , but with an "attached 1 ' 

n n 

velocity of Z . 



The combination of vortex sheet and the concentrated vortex in 
a fluid medium has to be force free since no external force is exerted 
on the system. 
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where the constants K are defined by: 
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NOTE: 



• » 

Since the M vortices are unattached, Z and Z are calculated 
from equation (29). 






Since there are more unknowns than there are equations, the system 
is still ill-defined. Four more equations can be gotten by imposing 
the Kutta condition at the point from which the vortex cores are fed. 
This is done by assuming that the fluid velocity relative to the cylin- 
der at these points is zero, i.e., Z ^ and are stagnation points. By 
solving equation (36) for -u^ + iv^ and for -u^ + iv^, then setting 
these equations equal to zero and rearranging in matrix form, the two 
new equations become; 
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The Ks are the same ones associated with equation (43). Though not in- 
cluded in equations (38) and (39) for the sake of simplicity, the above 
equations include the extra summation terms, and K^, necessary for 

the calculation if more than two vortex cores are present. To give the 
velocity equations (38) and (39) the same versatility, all that is neces- 
sary is to multiply these terms by i/llf* adding the first one to right 
hand side of (38) and the second one to (39). This has been done in 
the computer program for this model. 

The next two equations are gotten by solving the first equation of 
(42) for r\ and the second one for [ ^ • Then the time derivatives, 

and / , are taken, assuming that Z_ and Z, are constants. Since the 

z lb 

procedure is long and involved, but quite straight forward, only the re- 
sults will be presented; see equations (43). These equations also in- 
clude the terms necessary for more than two vortex cores. There is also 

* 

a provision made in these equations for U to be taken into considera- 
tion. This quantity was set to zero for the study undertaken in this 
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thesis, but could be utilized in future work. Note that the equations 

* 

have been non-dimensional ized . This means that U is also a non- 

oc 

dimensi ona 1 ized quantity. 

0 0 

This now leaves three six unknowns, , Z^, Z^, Z^* Z^, and Z^ 
in the system of equations without a corresponding number of defining 
equations. The latter two unknowns can be easily handled by forming 
the conjugates of both sides of (40) and (41) to yield. 



z, +(z,-z T )— = - 



and 
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(44) 



21 — Z 



as mentioned earlier, Z, and Z_ are considered constants. This is an 
assumption based on experimental observations at higher laminar 
Reynolds numbers. The origins of these feeding sheets are in realitv 

functions of time, but their motions are so small that the assumption 

• • 

seems justified. Besides, to consider Z, and Z anything but zero, 

b T 

would have made the derivation of (43) at least another order of magni- 
tude more tedious. This may, however, prove to be a necessary future 
sophistication of the model. As a result, these two locations are 
taken from experimental data to be about 50 degrees either side of the 
cylinder's downstream axis; i.e.,$ T = 50 and = -50 degrees. 

No mention thus far has been made concerning the propagation as- 
pects of this solution, however since time dependent solutions generally 
imply that starting and ending conditions are necessary, it is well to 
discuss these now with respect to Z^ and Z With the assignment of 
any values to these two parameters, there have been enough equations 
derived and assumptions made to have the proper balance between unknowns 
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and defining equations to describe the system at any given instant. Not 
just any value of and is desirable however, since this will not 
give even an approximate picture of core growth and drag buildup. Bryson 
[29], through certain mathematical considerations, found that the path 
line of the vortex cores, as they grow and move away from the cylinder, 
is unique in the sense that for any given feeding sheet origin, Z or Z^ , 
there can be only one path, per vortex, followed as time progresses. 
Without going into the details of his work, he also found that the feed- 
ing points, like the Foppl curve points,^ are mathematical equilibrium 
positions for symmetric vortex pairs. On the basis of these findings, 
it is also determined that cores can leave the cylinder's surface only 
in a direction 30 degrees downstream from the tangent to the cylinder at 
the sheet origin. Bryson arbitrarily selected the starting values for 
Z^ and Z ? as a distance of 0.03c away from the origins along these re- 
spective preferred lines of departure. This practice was continued in 
this study even though the model can handle asymmetry and more than two 
vortex cores. It was felt that these are good assumptions for an in- 
itial study of asymmetry. 

Except for determining when and how the vortex sheets are cut, there 
are no more assumptions to be made in completing the solution. The pro- 
cedure for solving the previously developed equations and the develop- 
ment of the (F + iF ) equation are not complicated and follow next. 

D L 



The Foppl curve is a symmetric curve about the downstream X-axis 
originating at the rear stagnation point, on which a vortex pair can be 
symmetrically located such that they will maintain their positions. This 
is much the condition that exists at Reynolds numbers below about 30. 

See [6, p. 223]. 
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Equations (40), (41), (43), and (44) can be rearranged and com- 
bined in a straight forward manner for the matrix equation (45). Both 
this matrix equation and equation (42) are solved using the partitioned 
matrix technique as given in Fox [39, p. 121]. Briefly, this involves 
breaking up the basic matrix equation into its real and imaginary parts 
and forming a new equation whose matrices are twice the order of the 
original. To describe the technique, suppose the original equation is; 



[a] [b] - {c} . 

where [A] is a square matrix of order j, while £b| , and^cj are column 
matrices with ] elements each. Since all the elements of A, B, and C 
are assumed to be complex, the matrices can be written as; 



w ■ 

W ■ { b r} + 1 [ b iJ • 

W ■ WW'i}- 

Substituting these expressions back into the original equation, results 

in 

<h] + i >J )( K} + i { b i]> ■ + 

Carrying out the indicated multiplication, and separating the real and 
the imaginary parts of both sides, yields 

NN - [NN • f c *j- 



and , 



K'l- 



b r] + l A J I B I j ■ { c i] 



These two equations are combined to form the new matrix equation: 

' N / \ 

H hi 



N -N 

IN N 
* — _ 









■; B ii f c i ; 



(46) 
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Since there are no imaginary numbers involved, this equation is solved 
by existing numerical methods, i.e., the inverse of the coefficient matrix 



premultiplies the constant matrix to yield the unknown values of 



In the computer program for this model, the inverse is found using a 
standard inversion subroutine with row pivoting. 

The propagation aspect of the solution is handled in the program 
through the use of a simple Runga-Kutta integration technique [40, p. 
232] by finding the subsequent vortex positions, after a time step, 

from equations (40) and (41). There are, of course, more sophisticated 
techniques, but in view of the parameters assumed, it is felt that the 
increased accuracy is not necessary at this point. In fact, there is 
not really that much difference in results if a straight Eulerian ap- 
proximation is used. 

The final relation to be developed is that for C + iC . This is 

U If 

done by completing the derivation of equation (30), or more specifically 
by finding the appropriate solution for equation (33). For this case, 
in which the cores are fed by a single feeding sheet originating at a 
point, Sarpkaya, [33], and [35], has previously derived the necessary 
equation, so the details will not be given here. The result, in non- 
dimensional form and adapted to this model is: 




c D +ic L =ziri 





(47) 
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The image velocities for the free and attached vortices are both calcu- 
lated, in a sense, using equation (37). However, where the results of 
equation (29) are used directly for the vortices, the results of solving 
matrix equation (45) must be used for the attached vortices. The modifi- 
cation to equation (37), used for this purpose is; for the image of 
vortex 1 , m 2? 



* 1 , =- 



and for that of vortex 2, 






.2 



(48) 



The computer program written for this model is listed in its 
entirety in the last portion of this thesis. It is written in FORTRAN 
IV as applied to the IBM 360/67 digital computer at the Naval Post- 
graduate School. In general, the calculations are carried out using 
sixteen digit numbers. The typical run, describing the impulsively 

started flow up to & 1 5 . 0 , requires less than three minutes of 

C 

computer execution time. Provided that the initial starting points, 
and are accurately defined, the solution appears to be stable and 

well behaved. The stability problem was previously discussed in con- 
nection with the vortex path line being unique. If the initial estimate 
for and Z^ are even slightly off, the solution technique attempts to 

return to the path, but as a consequence there may be some abnormally 

r l p-. ^ • 

^ ^ i* anc * ^2* T ° increase t ^ ie 

stability, smaller time steps are taken at the start and just after a 
vortex has been shed and a new vortex starts to grow. The normal time 
step is 0.05, while this reduced time step is 0.02. The program will 
accept as many as twenty vortices before it terminates execution. 
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Since anyone desiring to utilize this program must become intimately 
familiar with its internal workings, it is felt that any further dis- 
cussion is unnecessary here. There is a brief, but more detailed des- 
cription in Appendix A. 

It was possible, through the use of this program, to generate a 
reasonably good drag curve using asymmetric considerations. To start 
with t if no shedding was allowed, it was found that by multiplying the 

value of r { , during each time step, from equation (42) by 1.05, prior 

0 

to its use in equation (45), that the value of f 1 never changed sign. 

It did however, peak at about the same £ that it had previously 

(no asymmetry), but instead of decreasing to zero and subsequently chang- 
ing sign; it decreased somewhat and then increased beyond bound. The 
drag curve followed somewhat the same pattern. By arbitrarily cutting 
the vortex sheet at the minimum f~^ values, the drag curve in Figure 13 
was generated. The computer listing for Model A, though meant to be 
modified for other uses, is in fact, the same program that produced 
this curve. It is realized that the amount of asymmetry imposed is much 
too large, but is an illustration of the possibilities. A good begin- 
ning point for future work, would be with von Karman's stability condi- 
tions or to investigate the works done by Roshko [19] and Bearman [27]. 

D. MODEL B 

Since this model is an attempt to approximate the vortex sheets 
generated by boundary layer separations, it is much more sensitive to 
the placement of the nascent vortices than is Model A to the location 
of the feeding sheet origins. Hence, the problem of mathematically 
describing vortex sheet formation is more obvious. This model has a 
distinct advantage over the other, in that it is not necessary to 
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incorporate any assumptions into the derivation of the basic equations 
of motion that take the place of matrix equation (45). The assumptions 
are present, however, in the determination of the strength and initial 
nascent vortex positions. Probably the most important single advantage 
of this model, is that in general, once the separation points and vortex 
strengths have been determined, the vortex sheets will roll up to form 
vortex cores, taking care of the vorticity transport and diffusion auto- 
matically. To discuss these aspects further will require the develop- 
ment of the applicable equations, so this follows next. 

Since the nascent vortex is assumed to gain its total strength 
during the time step in which it is born, equation (29) can be used to 
determine its convective velocity. It is repeated here for convenience. 



-11+ A 



' V -u/£-,)+±[Z£-AJLl 



(29) Repeat 



The vortex convective velocity, as indicated in the initial discussion 

of equation (29), is calculated by setting Z = Z and ignoring the term 

n 

in the first summation that has the denominator equal to zero. The 
velocity of the image vortex can still be calculated from 



n, , ‘w -*7 l + +Vn. 

^ i U7L~ f ~' C r ZT2 



(37) Repeat 



■7i 



The drag and lift coefficient^ are again calculated from equation (30), 
however, this time, some liberties can be taken with its use. Again, 
since the nascent vortex is born during one time step, a part of the 
effect of the unsteady portion of (30), more specifically parts of equa- 
tion (33), can be assumed negligible. In fact, a study of the results 
of [41] indicates that the algebraic sum of the values of the terms 
concerned with feeding and growth of the vortices actually only account 
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1 2 . ’-** 



Configuration of Model B Showing tne 
Birth Areas and the Rotation Direction 
of the Vortices Born Therein. 




for about three percent of the drag. The two terras being discussed are 
equivalent to the terms inside the first and second sets of brackets on 
the right hand side of equation (47). Since in this case M * N, the rest 



of the terms become applicable to this model. Therefore, 




( 49 ) 



The system of basic equations is now complete, and like Model A, 
it has been programmed to allow asymmetry. The programming of equation 
(49) does not include the U_^ terra since during the study conducted 
for this thesis, only a uniform, steady cross-flow of * 1 was 

considered. It is simple enough to include this term, but then the 
deleted terras may become more important such that some kind of deter- 
mination might have to be made for them. It would be hoped that by 
the time these terms would be needed, the uniform flow problem would 
be solved and its results would be used to make this determination. 

This program is also written in FORTRAN IV, but unlike that of 
Model A, its calculation is carried out in only seven digits. The 
program, as written, will allow as many as one hundred nascent vortices 
to be born in any of the four birth areas during any given run. A run, 
which goes through the complete 100 time steps allowed, requires approxi- 
mately 30 minutes of computer execution time. Since the propagation 
aspects of the model are handled by Eulerian approximations, this could 
have serious error implications. However, since an approximate solution 
is all that is presently desired, it is felt that these errors are ac- 
ceptable. The large execution time is due primarily to the large number 
of vortices involved in the evaluation of (29), particularly at high 3st. 
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This must be done for each vortex, during each time step. The model also 
has the option of considering only the primary birth areas, 1 and 2, 
rather than all four. 

Before discussing some of the output from this model, the general 
method of determining vortex strengths must be considered. It has been 
shown by Fage and Johansen [42] that in the absence of back-flow and in 
a steady flow situation, the amount of vorticity leaving a separation 
zone can be accurately approximated by 



where V is the fluid velocity just outside the boundary layer at the 
point of separation. Thus, making the time steps, /\t, small enough, 
the nascent vortex strength can be given by 



The key problem now is to decide where and how to measure V and where 
the nascent vortex of this strength is placed. The four birth area 
model purports to take care of back-flow through the use of the second- 
ary birth areas, 3 and 4. Now since the maximum velocities on or near 
the cylinder's surface in these regions are very small as compared to 
the velocities in the primary regions, the effect of making some as- 
sumptions here is not too critical. The assumption made in most of the 
computer runs was that the maximum back-flow velocity on the cylinder's 
surface in that birth area was used in equation (51). From experimental 
studies, these back-flow vortices appear to separate at points located 




) 



(50) 




(51) 
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about + 55 degrees from the positive x-axis on the cylinder* s surface; 
hence, these are the assumed angles for the model. Instead of placing 
them on the surface, however, the model was generally programmed to 
place them at about 0.01c. This still leaves the problem of placing 
the primary nascent vortices; to the author *s knowledge, no truely 
satisfactory method does exist. The drag curves, and flow patterns 
shown in the figures of this subsection were all generated using arbi- 
trarily selected placement locations. The maximum velocity found on 
cylinder* s surface in these birth areas was, however, used to compute 
the strength from (51). For all the birth areas, the sign given to the 
resultant value of strength is as depicted in Figure 14. This sign con- 
vention comes from boundary layer considerations and the direction of 
the maximum velocity in the respective birth area. 

In order to prevent the extremely high convective velocities that 
can occur if and when two vortices get too close together, the model 
provides for coalescence and cancellation of the vortices depending on 
their comparative strengths. The vortices are combined if the alge- 
braic sum is not zero, in which case they are cancelled, by creating a 
new vortex having that sum for its strength and located at the "center 
of gravity'* of the two vortices. With a great many vortices this pro- 
cedure can consume much computer execution time. 

As in the case of Model A, any future user of this program would 
of necessity have to have an intimate knowledge of its inner workings. 
Appendix A contains a brief description of the parameters and the major 
sections of the program. 

This model was used to test a great many ideas using a great many 
runs. Figures 15 through 19 are presented to illustrate the kind and 
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type of output the model will generate* The program included with this 
thesis generated these figures. The primary vortices are injected into 
the stream at points whose angles from the positive x-axis are + 95 
degrees, and are located 0.045c from the cylinder's surface. 
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IV. DISCUSSION OF RESULTS AND CONCLUSIONS 



The primary result of this study is that reasonably representative 
drag coefficient curves can be generated using the two mathematical 
models developed. This is particularly true in the case of Model B t 
where the curves can be generated for both symmetric and asymmetric 
wake development. Figures 15 through 19 show the drag curve and wake 
development for an asymmetric wake case, while Figures 20 through 24 
show the symmetric wake case. The assumptions involved only the place- 
ment locations of the nascent vortices and the coalescence distance used 
to prevent abnormally high velocities if two nascent vortices come too 
close together. The symmetric wake case was generated by making the same 
assumptions for placement location as in the asymmetric case, however, 
the coalescence distance was increased. Where this distance was 0.06c 
for asymmetry, in the symmetric case it was 0.08c. These values were 
arbitrarily selected. 

In the case of Model B, the curve in Figure 15 was generated by 
assuming asymmetric generation of vorticity and assuming the vortex 



minimum. The amount of asymmetry and the time of sheet cutting is 
loosely based on the results of Sarpkaya’s experimental work [24]. 

Aside from the above results, numerous computer runs were made in 
an effort to develop a method by which the separation characteristics 
could be analytically described. Though these attempts were not en- 
tirely successful, some information was gained about the character of 
the vortex sheet origin and the requirements necessary to produce a 
reasonably realistic drag curve and wake structure. Of the two models, 



sheet is cut, freeing the core, when 




reached its first post-peak 
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Model B was the most illuminating since the concept of a string of small 
strength vortices is more nearly like nature than is the concept upon 
which Model A is based. As a result. Model B was utilized more exten- 
sively in this study. It is for this reason then, that the majority of 
the remaining remarks in this section will be concerned with results from 
this model. 

Of the various factors that go into the final description of the 
wake and computation of a drag curve, it was found that the location at 
which the nascent vortices are placed in Model B is by far the most 
critical parameter involved. Of the two items which determine the loca- 
tion of the vortex, angle and distance from the cylinder’s surface, it 
appears that the latter is more important. It was observed that there 
was a tendency for the primary vortices to "hug" the surface (Figure 25), 
all the way to the rear stagnation point, if their initial placement 
distance was too small. "Too small" in this case appears to mean less 
than about 0.03c. In Figure 25, the nascent vortices were placed at 
0.01c away from the cylinder’s surface, at angles of + 95 degrees. In- 
creasing or decreasing irj did not seem to have much effect on this 
tendency. Changing the secondary vortex birth location, for a fixed 
primary location, did not prevent this "hugging," nor did it appreci- 
ably change the roll-up tendency of the primary vortices if they were 
being shed reasonably well. Though the resulting wake formation in 
Figure 25 is not that desired for Reynolds numbers less than about 2 x 10"\ 
there is a striking resemblance to real flow patterns evident in the 
super critical Reynolds number region. By making approximately 50 com- 
puter runs in which only the birth location of the primary vortices was 
varied, it was found that a distance of about 0.045c and an angle of 
about + 95 degrees from the rear stagnation point, produced the most 
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reasonable flow patterns and drag curve; Figures 15 through 19. This 
shedding angle is comparable with the observed angle of separation of 
about 98 degrees (measured in the same manner). The distance from the 
cylinder's surface however, is about twice the value predicted by the 
Blasius series solution of the boundary layer equations. It is suggested 
therefore, that any future work in this area might well begin with a 
comprehensive study of the boundary layer transition region. A con- 
sideration of these regions above flat plates could be the starting 
point since the added problem of body curvature is not present. 

It does not seem possible, after studying the works of many in- 
vestigators, that Thompson [13] could have produced asymmetry with 
equations that are in themselves symmetric, were it not for the very 
real problem of computer round-off error. The fact that Figure 15 ex- 
hibits asymmetry would be of deep concern had Abernathy and Kronauer 
[10] not found that strings of small strength vortices would in fact roll 
up in an asymmetric manner. In the case of Model B, the predominate 
errors are due to assuming Eulerian propagation approximations and not 
from round-off. The Eulerian approximations are, of course, symmetric, 
so it is assumed that the asymmetry is in fact due to the interactions 
Abernathy and Kronauer observed. Interestingly enough, when using Model 
B asymmetry was not always produced, which leads to another conclusion; 
there is a critical combination of separation points and spacing dis- 
tance between sheets to cause the necessary interaction that forms the 
Karman street. On the other hand, reasonable drag coefficient curves 
could be generated without any asymmetry present. The overall conclusion 
is that the drag curve is predominatly influenced by the rate of vorti- 
city creation and its subsequent diffusion into the wake. This then, 
points right back to the over-riding problem of boundary layer separation 
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Since Model A feeds large vortex cores rather than creating vortex 
sheets, it would appear this model is better suited for describing the 
flow at Reynolds numbers less than about 50. By increasing the angles 
of the sheet origin, the drag curve peak can be increased to almost any 
desired value. In order to make the model applicable to the higher 
Reynolds number range, it will be necessary to find some scientifically 
based method of shedding these cores. The stability criterion of von 
Karman and others may provide the answers. 
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APPENDIX A 



A. DESCRIPTION OF COMPUTER PROGRAM FOR MODEL A 

The complete program for this model includes a primary or main pro- 
gram with seven subroutines that accompany it. There is one other sub- 
routine that is called that is not part of the listing included with 
this thesis. The subroutine is PLOTP and is called from the resident 
subroutine library of the Naval Postgraduate School Computing Center. 

The subroutine is used to produce printer plots, such as Figures 15 
through 19, instead of using the CALCOMP continuous line plotters as 
output, so it is not critical to the actual program calculations. The 
curve plotting subroutines or procedures of other systems could be 
easily substituted for this one. 

At the beginning of the main program and each subroutine is a short 
series of comment cards that indicate the purpose of the program and 
generally the parameters used therein. In general, the same notation, 
in FORTRAN format, is used in the program as was used to develop the equa- 
tions in Section III(c); for example in equation (43) is coded as K1 , 
or is coded as Z1B. Each computational section of the main program 
is preceded by comment cards indicating what computations that section 
is performing. 

Since there is no on-line function for extracting the imaginary 
part of a double precision complex number in the IBM version of FORTRAN 
IV, EQUIVALENCE statements had to be used to recover these values. 

It is assumed that the reader already has a working knowledge of 
FORTRAN coding techniques so that a definition list of the various para- 
meter names used in the program will be sufficient guide for future use. 
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The need for understanding the program is brought about by the fact that 
in order to study vortex sheet cutting procedures, the logic of the 
coding will have to be changed each time a different procedure is used* 
The section of this thesis entitled Computer Output, contains several 
pages of output that are to be used as trial run checks for a new pro- 
gram deck punched from the program listing* 

The type specification and precision of the following parameters is 
clearly indicated in the first few statements of the program, so there 
is no need for repeating it here* All arrays are stored in vector form, 
i.e., the values of the array elements are stored one column following 
another, going from left to right in the array. 

If a parameter is not found in the following list (equation para- 
meters described above are not included) but found in the program, then 
it is a dummy parameter used for intermediate operations only and there- 
fore has no special significance as compared to those listed. 



B. DEFINITION OF PARAMETERS 



AA(144) Vector containing the partitioned matrix elements. 



AC(12) Vector containing the partitioned elements of the right 

hand side of equation (45). First six elements contain 
the real parts, the second six the imaginary parts. 

AG(12) Vector containing the partitioned elements of the un- 

knowns in equation (45). 



CDCL(400) 

CV(20) 

DGl (401) 
DG2 (401) 



Vector containing the lift and drag coefficient values 
for each time step. 

Vector containing the free vortex velocities; corresponds 
to ZF (20) and GF(20). 

Contains for each time step. 

» 

Contains [~^ for each time step. 
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DG1MAX 
DG2MAX ^ 
DG1MIN 
DG2MIN J 


Contains the values of the maximum and minimum and [~^ 

for a run; used to determine when to cut the vortex 
sheets . 


DT 


Contains value of Zi t ; time step size. 


DUDT 


0 


EG 


-12 

A tolerance value, usually 1.0x10 


EGS 


EG squared. 


GF ( 20) 


Vector containing the values of the strengths ( ) 

of the free vortex cores; corresponds to the position 
vector ZF (20) . 


GG1 


Value of during a time step. 


GG2 


Value of P 2 during a time step. 


IT 


Counter value for the pass number through the program. 

1 ^ IT nit 


INNB 


If bottom vortex is shed during a time step, this para- 
meter equals 0, if vortex is not shed, it equals 1. 

It is used to determine whether a new vortex should be 
started from bottom sheet. 


INNT 


Has same indication as INNB, but applies to the top 
vortex and sheet. 


NCHDT 


Indicates the pass at which the value of DT is changed 
to prevent instabilities in the solution. 


NFV 


Counter used to indicate the number of free vortices 
existing. NFV ^ 20 


NFVT 


Contains the maximum value of NFV during any run. If 
NFVT >20, the DIMENSION statement must be changed as 
indicated in the program. 


NIT 


Contains the maximum value of IT during any run. If 
NIT > 400, the DIMENSION statement must be changed as 
indicated in the program. 


NKER(2) 


Vector containing the matrix singularity determination 
from subroutine INVERT. NKER(l) applies to equation 
(42) while NKER(2) applies to equation (45). 


NPL 


Indicates number of passes through program between wake 
development plots . 
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NPRT 


Indicates number of passes through program between 
status printouts. 


PI 


7r 


PI2 


7T 2 


T(401) 


Contains value of total elapsed for each time step 

c 


ZA(36) 


Vector containing the values of the elements in the co- 
efficient matrix of equation (45). 


ZAA(4) 


Same as ZA(36), but applies to equation (42). 


ZC(6) 


Vector containing the elements of the column matrix on 
the right hand side of equation (45) . 


ZCC(2) 


Same as ZC(6), but applies to equation (42). 


ZCP(IOO) 


Contains the circle circumference points and scaling 
points used for plotting purposes. 


ZF(20) 


Contains the free vortex positions during a time step; 
corresponds to the vector GF(20) on a one for one basis. 


ZG(6) 


Contains the values of the unknowns in equation (45). 


ZGG(2) 


Same as ZG(6), but applies to equation (42). 
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APPENDIX B 



A. DESCRIPTION OF COMPUTER PROGRAM FOR MODEL B 

This program is much shorter from a coding standpoint than is Model 

A. There is only one subroutine besides the plotting routine, PLOTP, 
needed for operation of the program, and it is included in the listing. 
The general comments previously made for Model A, apply equally well 
here. The only major difference, besides the concepts upon which the 
two are based, is that this program is written using single precision 
arithmetic rather than double precision as in the case of Model A. 

B, DEFINITION OF PARAMETERS 



AVN(4) 


Values, in degrees, of the nascent vortex placement 
angles actually used in the calculations. AVN(l) 
applies to birth area 1; AVN(2) applies to birth area 
2, etc. 


CDCL(IOO) 


Same as Model A. 


CLOCKR ~) 
CLOCKS V 
CLOCKT J 


These parameters are used to determine the actual 
elapsed computer time for each pass and run; 
can be eliminated from the program without hurting 
its purpose. 


CW (400) 


Contains the value velocities for each vortex during 
a time step. CW(l)to CW(100) Birth Area 1 

CVV(lOl) to CVV(200) Birth Area 2 
CW(201) to CW (300) Birth Area 3 
CW(301) to CW(400) Birth Area 4 


DTP 


Same as DT in Model A 


EG 


Same as Model A 


EGS 


Same as Model A 


GK(400) 


Contains the value of the vortex strengths during a 
time step; corresponds to CW(400) on a one for one 
basis . 


IT 


Same as Model A 
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ITIME(O) 


An IBM on-line function used to make real time 
available to a FORTRAN IV program; can be eliminated 
without hurting the purpose of the program. 


MPA 


The value of nascent vortex placement distance as a 
function of RE as generated by formulas from Sarpkaya's 
paper [35]. Used for comparative purposes. 


MP(4) 


Contains the actual nascent vortex placement distance 
from the cylinder's surface used in the calculations. 
MP(1) corresponds to AVN(l) or NVM(l) ; MP(2) corres- 
ponds to AVN(2) or NVM(2) ; etc. 


MS 


Radius of the arc along which the search for maximum 
velocities is made; usually equal to 1.0. 


MZD 


Minimum distance two vortices can be from each other 
without being coalesced or cancelled. 


NB 


Total number of birth areas; usually equal to 4, but 
if only the primary nascent vortices are being studied, 
then NB=2 . 


NCAN 


Counter indicating number of vortices cancelled during 
a run. 


NCL 


Counter indicating number of vortices coalesced during 
a run. 


NIT 


Same as Model A. 


NP1 


Indicates number of passes, times ten, between wake 
development plots. 


NV 


Counter indicating number of primary vortices generated 
in birth area 1 during a run. It is also equal to the 
number generated in birth area 2. 


NVM(4) 


Contains the values, in degrees, of the maximum velocity 
locations resulting from a search for same in the birth 
areas. NVM(l) corresponds toAVN(l); etc. 


PI 


Ir 


RE 


Reynolds number; used in conjunction with MPA. 


TP(lOl) 


Same as T(401) in Model A. 


VELM(4) 


Contains the values of maximum velocity found during the 
search for same in the birth areas. VE12i(l) corresponds 
to NVM(l) ; etc . 


ZCC(IOO) 


Same as ZCP(IOO) in Model A. 
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ZCFP(51) 


Contains the cylinder* s circumferential locations from 
135 to 85 degrees (measured from the rear stagnation 
point) . 


ZCRP(61) 


Contains the cylinder* s circumferential locations from 
0 to 60 degrees (measured from the rear stagnation point) 


ZF (400) 


Contains the locations of the vortices during a time step 
Corresponds to CW(400) and GK(400) . 
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COMPUTER OUTPUT, MODEL 
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